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Abstract 

To learn about basic aspects of nano-scale spherical molecular shells 
during their formation, spherically curved two-dimensional N -particle 
Lennard- Jones systems are simulated, studying curvature evolution paths 
at zero-temperature. For many N -values (N < 800) equilibrium configu- 
rations are traced as a function of the curvature radius R. Sharp jumps 
for tiny changes in R between trajectories with major differences in topo- 
logical structure correspond to avalanche-like transitions. For a typical 
case, N = 25, equilibrium configurations fall on smooth trajectories in 
state space which can be traced in the E — R plane. The trajectories 
show-up with local energy minima, from which growth in N at steady 
curvature can develop. 

1 Introduction 

In vitro self-organization in aqueous solution of nanoscale spherical shells, 
like various types of nano- vesicles ^ 121 El EI) and viral capsids is a ther- 
modynamic process driven by overall free energy minimization. The result- 
ing shell has inherent tendencies toward a) crystalline structure [HI El j b) 
global polyhedral svmmetrv [HI ITUl [TT] . and c) discrete sizes ^^1^1, thus 
discrete curvatures. 

A compelling question is whether during thermodynamic growth a tran- 
sient uncompleted (open) shell of given N prefers intrinsically to adopt a 
specific equilibrium radius of curvature (i?), and whether this radius - on 
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the basis of molecular packing order already - may be approximately stable 
during (part of) that growth. The occurrence of local minima in the mean 
(per particle) potential energy (E) of an emerging shell as a function of all 
internal degrees of freedom - including R - along the evolution path in state 
space, could thus be significant for steady curvature during growth in A^. 

In this paper we approach the problem in a much simplified model that 
allows systematic generic studies. The sphericity of the shell, in reality due 
to intrinsic SD-properties of the optimally closely packed molecular subunits, 
is built-in as a global geometrical constraint. We perform computer experi- 
ments with freely relaxing - possibly open - spherically curved zero-thickness 
monolayers of identical molecules, studying structural properties in relation 
with N and R. Our approach is most natural for monolayer nanowheel vesi- 
cles jl], but it is as well relevant for tiny bilayer vesicles, where the laterally 
most densely packed pseudo-crystalline inner lipid headgroup-sublayer can 
act as monolayer 'backbone'. 

As model we use two-dimensional (2D) Lennard-Jones (LJ) A^-particle 
systems on a spherical surface with flexible radius at zero temperature jHl 
in] . The LJ potential Vij = r~^'^^^ — 2rj^^ between two particles i and j 
with Euclidean distance rjj energetically favors close regular packing with 
essentially unit distance between neighboring particles [H]. The LJ-form in 
the constrained system acts as an effective interaction mimicking the real 
complex of interactions. It allows a comprehensive systematic exploration 
while keeping salient features of real systems. The present work arises from 
a series of computational studies 1^ on 2D spherical crystallization in LJ- 
systems over a broad range of A^, involving thermodynamic behavior and 
zero-temperature global energy minimization. The LJ-systems follow local- 
equilibrium paths in state space, realistically allowing for local minimum 
'hang-ups' in evolving configurations. 

When many LJ particles are randomly spread over a flat surface they ag- 
gregate into an approximately homogeneous configuration, a major fraction 
of the particles being trapped inside the bulk (interior). Edge particles have 
higher energy than bulk particles, giving rise to edge tension, the 2D equiv- 
alent of surface tension. Minimizing edge energy, flat aggregates become 
approximately circular patches. Minimizing the overall potential energy the 
bulk becomes an essentially regular hexagonal lattice. Allowing for spherical 
curvature, the 2D system can further decrease the energy by reducing the 
edge length and by a rising attraction from LJ-tails of remote particles (the 
LJ-forces acting in 3D). The energy gain by curvature, however, balances 
against increasing strain energy of the bulk because of less favorable packing. 

Our question becomes whether - thanks to the interplay of these in- 
dependent non-linear effects - equilibrium radii Req can be found where 
relaxation after any small change in curvature raises and freely variable 
R leads the system back to the same equilibrium radius. If such local energy 
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minima in open configurations are thermodynamically significant, tliey can 
stabilize transient states along a path of growing A^. This is supported 
by a L.T-studv |15j. showing that closed global minimum energy A^-particle 
configurations (covering the whole sphere) strongly correlate with specific 
open local minimum configurations of lower A^. 



2 Methods 

In our experiments, each time the radius R is changed by a small step, the 
particle system is relaxed by minimizing - at the new R value - the mean 
energy, which implies reaching the nearest equilibrium configuration. The 
system is thus evolved in curvature while staying in equilibrium with chang- 
ing R. For relaxation (equilibration of forces and energy minimization), 
a steepest descent (SD) algorithm is applied while for aggregation a 
Metropolis Monte Carlo simulated annealing (SA) optimization is used. 
The latter method enhances the probability that the system ends up in a 
global energy minimum (GEM) jTH] rather than in a local secondary energy 
minimum. The 2D-topology of a configuration is defined by the Voronoi 
nearest neighbors method ) giving each particle a coordination number 
(CN), which is 6 everywhere for a flat (hexagonal) GEM configuration. The 
value CA^ — 6 is denoted as the disclination charge - short d-charge - of a par- 
ticle. In a 2D topological structure any single built-in defect - disclination 
or dislocation (tightly bound pair of disclinations with opposite d-charge) 
- can be displaced but not removed, except by its annihilation as part of 
a set of converging complementary defects, or by moving it all the way to 
the edge. During the transformation of a flat GEM lattice to a closed shell 
a net total d-charge of -12 must be incorporated in the full Voronoi lattice. 
In addition, dislocations have the function of lowering the strain energy in 
total, by distributing it locally more evenly . 

How much the configurations change during relaxation after a step AR, 
can be expressed as the mean Euclidean distance, f, traveled between the 
associated sets of coordinates xjj and :x-r+ar in configuration space (mean 
taken per particle and per percent change in curvature): 




3 Computer experiments 



From flat to spherical 
In a first series of experiment we explore gross changes in topological 
and geometrical structure and the corresponding energy with monotonously 
decreasing R. For hundreds of runs with A^ < 800 an initial, circularlike 
aggregate is prepared from a flat regular hexagonal lattice with unit spacing 
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(GEM for infinite in flat 2D). In the experiments R is decreased in 1% 
steps, each time the configuration being projected onto the new sphere and 
then relaxed using the SD method. The decrease in R is continued until 
the system is compressed considerably. A system of particular follows a 
'standard' (for that A^) evolution path through state space and through the 
E — i?-plane. 

Example, N = 25 

Secondly, curved lattices unbiased by any initial configuration and path 
history, are simulated while starting in SA mode at high temperature (T = 
10) in a random configuration and then aggregate by cooling down the sys- 
tem in 5 % AT/T steps. SD is applied as final tuning. This study is done 
for a typical 'unmagic' N value, = 25. In 1300 runs radii are randomly 
chosen between R = 1.2> and R = 2. 

The equilibrium points va. E — R space align over a range of i?-values 
along distinct smooth trajectories: lines which correspond to continuous sets 
of topologically and geometrically closely similar equilibrium configurations. 

Tracing up and down in curvature 
The central question of the current paper, whether trajectories of the 
system lining-up closely related configurations, can provide a stable system 
against freely variable -R, is addressed in a third type of experiment. Starting 
from particular open A^ = 25 configurations obtained in the second exper- 
iment, trajectories are traced (with relaxation) step- by-step in R in both 
directions, and the structure is investigated. 

4 Results and discussion 

Energy trends 

In the first study a general energy trend is found for all A^- values, as 
illustrated by three typical examples: for N=6, 50 and 350 in Fig.^ (left 
column). As expected, the most prominent feature is a deep global mini- 
mum along the followed E — R path. Having started from a flat regular 
configuration while systematically decreasing R, the particles cover in close 
packing an increasing part of the sphere until at the closure radius Rc the 
lowest E along the E - R path is reached. At this stage any uncovered area 
- and thus any edge - has disappeared. We note that except for details, the 
structure is not biased significantly by the initial configuration. 

The energy difference, AEc, between the — > oo asymptotic value and 
the deep minimum represents the energy gain of closure. Additional struc- 
ture is visible: secondary, local, minima and sharp downward jumps at 
R > Rc. 

Toward smaller radii, R < Rc, the system comes under external pres- 
sure, and the density increases. The LJ repulsive core between neighboring 
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Figure 1: Data from the first series of experiments: from top to bottom for = 6, 
50 and 350 the mean potential energy (E) (at left), and the mean Euclidean distance 
{M ED, f) traveled in configuration space during equilibration after 1% downward steps 
in curvature radius {R). The system was initiated in a flat regular configuration. 



particles then increasingly dominates the energy, the energy rises steeply, 
independent of the detailed configuration. 

The energy of a configuration can be approximated as E = (-/VbuikCbuik + 
-^edgcecdgc)/-^, where N^^u, + N^dge = N, ebuik = -3.382 being the lowest 
possible value corresponding to an infinite flat regular hexagonal lattice while 
eedgc is typically half that value. Due to the impossibility of fully regular 
hexagonal packing ebuik is higher on a curved surface. The unfavorable 
edge energy forms the main driving force for curvature and shell closure, 
although - due to topological rearrangement - normally a barrier between 
fiat and curved states stands in the way. Consequently, for the larger A^- 
values, starting from the flat GEM configuration, indeed dE/dR is negative, 
meaning that at sufficiently low T such a flat LJ-system is locally stable 
against curvature. An exception is the extreme case, N = 6, where the 
initially flat LJ-system can immediately gain energy by curving. A tiny 
jump in E (at R ^ 1.36) goes with a discontinuity in the dE/dR slope. 
Here already we observe a secondary minimum along the standard E — R 
path. For iV = 50 below i? = 4, local energy minima and jumps are visible. 
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Figure 2: Mean energy {E) versus curvature radius (R) for N = 25. From top to 
bottom: a) Second experiment: energy minima at random fixed values of R, found 
in unbiased searches, each with random initial configuration, applying SA and SD; b) 
Third experiment: example of tracing through E — R space by 1% steps in R using SD. 
Arrows indicate the R step direction. The inset shows a clear-cut case of hysteresis; c) 
Parts of trajectories, exhibiting local minima as a function of R. 



At Rc ~ 1.87 the system closes over the sphere. For = 350 below R = 13, 
many minima and many jumps show-up. At Rc ^ 5 the system closes. In 
the second experiment {N = 25), in Fig. [2^, the independently obtained E- 
values for fixed, randomly chosen, i2-values align along smooth lines in the 
E — R plane. 

Euclidean distance 
Energy jumps (Fig.l at left) go with structural transitions visible as 
spikes in the mean Euclidean distance per particle per percent, f{R) (Fig.l 
at right). At -R < Rc-, the particles hardly move anymore over the sphere, 
all flexibility being lost due to the external pressure counterbalanced by the 
repulsive core. The smooth background in r[R) reflects small adjustments 
of an essentially stable configuration after a curvature step. It can be fully 
suppressed by reducing the step size, which turns out to leave the transitions 
essentially unaffected; the sharpness of the transitions remains within a step 
size of even an order of magnitude smaller than applied in the current data. 
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Closure radius 

The radius Rc at which the system closes can be predicted by requiring 
the packing density on the spherical surface to be (almost) the same as for 

a flat lattice, Rc{N) = y^£:Cos(f ). For iV=6, 25, 50 and 350 it follows that 
i?c(-^)=0.64, 1.3, 1.9 en 4.9 respectively. These values agree with the data 
(Fig. 1 left) except for small N where the packing is far less optimal than 
for a flat GEM lattice. 

Structural transitions and topological defects 
Energy jumps between different trajectories involve a major global rear- 
rangement. 

With successively applied small steps in curvature followed by relaxation, 
the system usually undergoes only minor local redistributions, while strain 
accumulates in parts of the lattice. At some curvature steps, however, an 
avalanche of sequential displacements over a major part of the system is 
triggered by local instabilities, while releasing much of the built-up strain. A 
small change in R may thus enforce a structural transition to a conflguration, 
in which both the strain and the 2D-topological (defect) structure have been 
drastically altered. An extensive systematic study of the topological defect 
structure as a function of up to 200 is available ^6^ . 

The transitions come along with incorporation or removal of defects as 
a function of curvature. Energy barriers are present between configurations 
with different topological structure. The transition to take place requires 
the system to climb the barrier to a threshold for starting an avalanche of 
particle moves. The threshold ii-value depends on the barrier side, causing 
hysteresis. 

Starting from an essentially circular patch with flat regular packing, the 
'unstable' radius Ru where - with increasing curvature - the flrst structural 
transition occurs, is modeled as S * 2nmax ai'csin(a) = sin(2nmax arcsin(a)) 
[5], where a = l/(2i?u), ?^max is essentially the largest completed hexagonal 
ring, and where S can be taken from the data for a single speciflc A'^-value 
(see Table^. 

Open configurations 
The points in Fig.|2^ at large R cover a broad range in energy, corre- 
sponding to a great variety of open conflgurations and edge-arrangements. 
As R decreases, the edge becomes smaller, and the variation decreases. Any 
flexibility essentially disappears below R^ due to the strong constraints for 
a closed spherical configuration. 

Secondary minima 
Fig. El) shows E — R trajectories from the third experiment, starting from 
a specific configuration (near the minimum at R ~ 1.5) from the second 
experiment. Here the system is traced up-and-down in R. Indeed, the 
trajectories in the E — R plane connect unbiased solutions, and - like in 
the first experiment - jumps and secondary minima occur. The secondary 
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1.37 


1.34 


50 


3 


3.95 


4.11 


350 
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12.3 


12.3 


500 


11 


15.6 


15.0 



Table 1: The radius i?u where the trajectory becomes unstable for increasing curvature 
starting from a flat regular distribution for iV=6, 50, 350 and 500, compared with 
estimates from the simple model with S calibrated at = 350. 

minima visible in Fig. 2c, are obtained starting from configurations of the 
second experiment (Fig. 2c), and tracing in R up and down until a jump 
occurs. The minima have a typical depth of 10% of the closure energy AE'c. 
At non-zero T the thermodynamic significance of such minima should be 
judged with respect to both kT and AE'c. 

Variable N 

The finding of distinct locally-stable open configurations for fixed N 
raises the expectation that any discrete curvature may remain, or change 
smoothly, during growth when particles are added. A study of this type ^3] 
indeed indicates that discrete locally-stable configuration trajectories exist 
as a function of both R and N ^ where the bulk packing remains essentially 
the same. 

In conclusion, for fixed- two-dimensional spherical Lennard-Jones sys- 
tems at zero temperature, the global energy minimum with decreasing cur- 
vature radius R is approached through sharp transitions with major rear- 
rangements. These transitions bring- in topological defects connected with 
curvature. The curvature range of events relevant for self-assembly, Ru - Rc 
is consistent with simple models. Apart from the closed (global) minimum 
energy configuration, secondary, local, minima show up at larger ii- values, 
with an open configuration. This phenomena - here shown for A = 25 
- occurs naturally as a consequence of optimal packing topologies. Dur- 
ing growth toward a complete shell such minima can capture 'threads' of 
steady curvature along with growth in A^, while staying in equilibrium. The 
present results and methods can help guide further generic studies of the 
self-organization in complex spherical molecular systems. 

The authors are grateful to D. Frenkel (FOM-AMOLF) for many fruitful 
suggestions, and to M. Livny (UW-Madison) for indispensable support on 
Condor High Throughput Computing. Part of this work has been funded 
by the FOM and NWO science organizations in The Netherlands. 
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